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1 INTRODUCTION 



The discovery of the first Einstein rings at both radio fae witt et alj 
[l988) and optical ( Warren etal. 1996) wavelengths brought with 
them the promise of a new age of precision galaxy mass measure- 
ments. Strong lensing is a powerful astrophysical tool with applica- 
tions ranging from measuring galaxy structural properties and evo- 
lution to cosmological parameter measurements using time delays 
(for Ho) or the statistical properties of lens samples. 

Some applications of lensing, such as determining Bubble's 
constant (Ho) from time delays of multiple-image quasars, rely on 
accurate knowledge of both the radial mass profile of the lens mass 
and the location of the lens centre. Quadruple images of lensed 
QSOs commo nly have a range of solutions fo r mass slope and 
ellipticity (e.g. IWambseanss & Paczvnskilll994l) which cannot be 
resolved without additional assumptions or information. Einstein 
rings can, in principle, help to solve this problem because they 



ABSTRACT 

We have developed a new software tool, Lens VIEW, for modelling resolved gravitational 
lens images. Based on the LensMEM algorithm, the software finds the best fitting lens mass 
model and source brightness distribution using a maximum entropy constraint. The method 
can be used with any point spread function or lens model. We review the algorithm and intro- 
duce some significant improvements. We also investigate and discuss issues associated with 
the statistical uncertainties of models and model parameters and the issues of source plane 
size and source pixel size. 

We test the software on simulated optical and radio data to evaluate how well lens models 
can be recovered and with what accuracy. For optical data, lens model parameters can typically 
be recovered with better than 1% accuracy and the degeneracy between mass ellipticity and 
power law is reduced. For radio data, we find that systematic errors associated with using 
processed radio maps, rather than the visibilities, are of similar magnitude to the random 
errors. Hence analysing radio data in image space is still useful and meaningful. 

The software is applied to the optical arc HST J15433H-5352 and the radio ring 
MG1549H-3047 using a simple elliptical isothermal lens model. For HST J15433H-5352, the 
Einstein radius is 0'.'525 ±0.015 which probably includes a substantial convergence con- 
tribution from a neighbouring galaxy. For MG1549H-3047, the model has Einstein radius 
I'.'IOS ±0.005 and core radius 0'.'16 ±0.03. The total mass enclosed in the critical radius 
is 7.06 X lO^^M© for our best model. 

We finish by discussing issues relating to modelling of resolved lens images for this 
method and some alternatives. 

Key words: gravitational lensing - galaxies: individual: MG1549H-3047 - galaxies: individ- 
ual: HST J15433H-5253 - methods: numerical. 



roughly trace the shape of the tangential critical line, hence the el- 
lipticity of the mass distribution is well determined. In addition, the 
many lines-of-sight provided by Einstein rings sample the gravita- 
tional potential much more thoroughly than QSO lenses. Hence, 
there is a great deal of information contained in Einstein ring im- 
ages. 
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Modelling resolved gravitationally lensed images, such as 
Einstein rings, is a much more complicated problem than mod- 
elling lensed QSO images. The lensed image of a resolved source 
is determined both by the properties of the lensing galaxy and 
the source brightness profile. Consequently, algorithms for mod- 
elling resolved images have gone through several stages of de- 
velopment. The Ring Cycle algorithm (Kavser & Schramm 198^ 
lKochaneketalJll989t) reconstructed a source based on the as- 
sumption that the data is a true representati on of the lensed 
surface brightness. LensCLEAN 0Cochanek & NaravanI Il992h 
uses a CLEAN style algorithm to build a non-parametric source 
from a processed radi o map. This was subsequently improved 
iEllithorpe et alJ Il99d) to work directly on radio visibilities 
rather than CLEANed radio maps. More recently, it was im- 
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i iroved ag ain i\Vuc knitzl20o'4l) and applied to the lens B02 18+357 
Wucknitz et alJ2004 to make an accurate Hq mea surement. Also 
recent ly, the 'semilinear' method was developed bv lWarren & Dvd 
i2003h who noted that for a fixed lens model, determining the 
source model that best describes the data can be posed as a lin- 
ear inversion problem and c an therefore be tackled with sta ndard 
matrix inversion techniques. iBrewer & Lewis! i2005l l2006h have 
investigated the lens modelling problem using genetic algorithms 
and Bayesian techniques, respectively. 

An alternate approach to modelling radio images was provided 
by the LensMEM algorithm ( Wallington et al. 1994 , 1996, here- 
after WNK94,WKN96). This approach used a pixelised source and 
the novel idea of a mapping matrix (relating pixels in the image 
and source planes) to generate a source model using an iterative 
approach. This algorithm is equally suitable for radio and optical 
data with a point-spread function (PSF) and it forms the basis of 
the software described in this paper. 

In ij2| the algorithm is reviewed and some significant 
impr ovements ar e desc ribed, in particular those proposed by 
Iskilling & Brvai] [19841 hereafter SB84). Also discussed is the is- 
sue of defining an acceptable goodness-of-fit with non-parametric 
(i.e. pixelised) sources for both radio and optical data. In J^lthe 
software is applied to test optical and radio data, showing that for a 
radio map the lens model parameters can be recovered quite accu- 
rately, despite the problems associated with analysing CLEANed 
images. In addition, it is shown that error estimates for the lens 
model parameters can be reliably made for optical and radio im- 
ages using the properties of the surface. In ij4] LENS VIEW is 
applied to the optical arc HST J15433-I-5352 and the radio ring 
MG1549-^3047. In some issues relating to the software, and 
modelling of resolved lenses in general, are discussed and com- 
pared to some of the other methods outlined here. 

Throughout this paper, we use a flat, A-dominated cosmology 
with = 0.3, Q.A = 0.7 and Ho=lQ km s"^ Mpc"\ 



2 HOW LENSVIEW WORKS 
2.1 A review of LensMEM 

In this section the basic ideas behind the LensMEM algorithm 
are reviewed and some improvements and additions found in 
LENSVIEW are introduced. 

Lensview is essentially an iterative method for calculating a 
source brightness distribution that best matches an observed lensed 
image for a given lens model. LENSVIEW requires a mechanism to 
allow the projection of a resolved source model into a lensed image 
while correctly preserving surface brightness. The iterative nature 
of LENSVIEW also requires the ability to 'reverse-project' from the 
image plane back to the source plane (to provide feedback for the 
source model), an operation which is more complicated than for- 
ward projection because of the many-to-one relationship between 
image pixels and some source pixels. 

Projection and reverse-projection are achieved by creating a 
'mapping matrix' which relates each pixel in the image plane to 
each pixel in the source plane via the lens model. Pixels in the im- 
age and source planes represent surface brightness. To generate the 
mapping matrix, image pixels (denoted hj) are projected back to 
the source plane and the area of overlap for each source pixel (de- 
noted /fcj) is calculated. The image pixels are divided into two tri- 
angles so that the area is calculated correctly if the pixels are folded 
or twisted by the projection process. This happens to pixels which 



overlap critical lines. The area overlap, or 'weight' (denoted Wijki) 
measures what fraction of source pixel fki is covered by the pro- 
jected image pixel lij. The weight matrix Wijki has many useful 
properties which are repeated here from WKN96. 
A source can be projected into an image using 



(1) 



and is typically convolved with a point spread function or "beam", 
Bpq, to generate the model image, Mij, that a telescope would see. 
The inverse magnification of an image pixel fj,~^ is 



'Eki "^'Jfc' 

image 



(2) 



where Aimagc is the area of an image pixel in units of source pixel 
area. The magnification in the source plane fiki is 



(3) 



where the results of the summation inside the parentheses can be 
performed in advance and stored. 

In LensMEM, the software proceeds in two nested cycles: 
an inner loop and an outer loop. The inner loop finds the source 
that best fits the data (subject to the MEM constraint) for a fixed 
lens model. The outer loop adjusts the lens model parameters and 
begins the inner loop again. 

Mathematically, the aim of the inner loop is to minimise the 
function^ J = C — XS, where 



(4) 



is the usual merit function evaluated between the data, Dij , and 
model image and 



(5) 



is the entropy in the source model. A is the default pixel value for 
the source, or "sky background" (see Section l231 for further discus- 
sion). 

For each inner loop iteration, the source is incrementally up- 
dated /^;+^ = + eAfki. In LensMEM, the conjugate gradient 
method was used so Afki = VCki — XVSki + Pki where VSki is 
the gradient of the entropy^, VCki is the gradient of the reverse- 
projected back to the source plane and Pki is a weighted sum of 
previous A/^fs as defined by the conjugate gradient technique. 
Generating VCki is a complicated process and the details are dis- 
cussed in Section |2!4| Once Afki has been calculated, it must be 
scaled so that J is minimised for the iteration. To determine the 
scaling factor, e, a trial image using the source fki + eAfki is gen- 
erated, forward projected and convolved with the PSF to calculate 
X^- Thus at each inner loop iteration, the problem is reduced to a 
one-dimensional line minimisation involving several forward pro- 
jections to find the scale factor"', e. 



^ Note the different definition of J compared to WKN96. 

/fci, A/fej, VCfej, \'Ski and Pki ai'e images the size of the source plane. 
e and A are scalars. 

This line minimisation procedure must also be iterative because the gen- 
eration of the trial source, fki + eAfki, can, in principle, generate negative 
source pixel values which must be corrected. If the entropy constraint is not 
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The inner loop stops when the minimum in J is reached. At 
this point, A is adjusted and the inner loop re-started. This pro- 
cess continues for up to 100 iterations. WKN96 describe the A- 
adjustment process as beginning with entropy-dominated source 
models (large A in our terminology) then decreasing A until a good 
fit is found, without providing specific details. 



2.2 Departures from LensMEM 

Since this is a computational method, we are interested in algo- 
rithm accuracy, stability and efficiency. In the computational con- 
text, there are some issues with the implementation of LensMEM 
as defined in WKN96. Firstly, there is the process of adjusting A. 
As described, the inner loop is inefficient — effectively, a 'mid- 
dle' loop has been introduced where A is adjusted when the inner 
loop is making little progress, but before the lens model is changed. 
At the heart of this issue is a lack of well defined stopping point, 
or "good fit" to the data. We elaborate on this in i)2.3l From an 
algorithmic perspective, we want to adjust the source and A simul- 
taneously on the way to minimise J which reduces the number of 
forward/reverse projections required. 

Secondly, VCfcj and VSki are vectors that are the result of 
operations on completely different spaces, so adding them to form 
A/fei can lead to problems such as negative source pixels. In ad- 
dition, the definition of the entropy S is WKN96 allows individual 
pixels in the source to dominate both S (hence J) and VSki, which 
affects both stability and accuracy and is undesirable. 

Lastly, the fixed number of iterations of the inner loop artifi- 
cially limits the ability of a lens model to reproduce the data. Our 
experiments showed that models with a broad PSF (e.g. radio data) 
generally converge more slowly than data with a narrow PSF (e.g. 
HST data). Furthermore, WKN96 state that "with a good model, 
the final can typically be far smaller than A'^pix". This highlights 
the importance of knowing the target for the model. The point of 
MEM is to find the most likely source given the data subject to 
reaching the target value. Data with appropriate errors should not 
generate a significantly lower than the target even for the cor- 
rect model. Poor lens models are never able to reproduce the data, 
so will converge to some x^ greater than the target regardless of 
A. Using the Skilling & Bryan method, good lens models converge 
naturally to the target x^ because the entropy weight A is adjusted 
as part of the inner loop. With LENSVIEW, a more flexible criteria 
for stopping the inner loop is adopted. If the model does not show a 
fractional improvement of ~ 10^"^ between iterations, or the target 
X^ has been reached, the inner loop is stopped. Of course, one may 
want to find the best possible source given a particular lens model, 
in which case the target x^ can simply be set to some impossibly 
low value. In that case, the entropy effectively becomes a constraint 
that source pixels cannot be negative. 

SB84 identified these issues (and others) as general features 
of this kind of problem and proposed an elegant algorithm to solve 
the problem. In LENSVIEW we dispense with the conjugate gradi- 
ent method, the A-adjustment process and the stopping criteria of 
LensMEM and adopt SB84's general purpose algorithm to solve 
the minimisation problem. In Section 12.51 the Skilling & Bryan 
method is briefly reviewed and its implementation in LENSVIEW 
is described. 



u.sed (i.e. A = 0), then source pixels can be negative and multiple forward 
projections are not necessary. 



2.3 Determining the target x^ 

Determining the target x^ is essential if meaningful comparisons 
between lens models are going to be made. In the general case, the 
number of degrees of freedom, v, in a modelling problem is given 
by 

V = A^data — NmodA (6) 

where A'data is the number of independent data points and A^modci 
is the number of model parameters. For large v, the target x^ is 

V ± \/2v for a 68 percent confidence interval. However, a pixelised 
source model and (potentially) broad PSF creates ambiguity in u 
for the lens modelling problem. 

Kochanek ( 1995) considered the issue of degrees of freedom 
for modelling of CLEANed radio data and concluded that 



^bcam 

where Atg.n is the area inside the tangential caustic for the lens 
model, j4beam is the area of the beam and m is the number of pa- 
rameters in the lens model. In this case, x? is evaluated over the 
entire multiply-imaged region (singly-imaged regions can be fit- 
ted perfectly with point sources and provide no constraints). This 
definition is interesting because it is independent of the number of 
pixels (or CLEAN components) in the source plane. It is most use- 
ful for radio rings where a large fraction of the multiply imaged 
region contains lensed flux. For images where the flux is confined 
to a small region (usually around the tangential caustic) this is not 
a good definition because the many pixels that contain no flux are 
still counted in x^, thus reducing the power of the statistic to dif- 
ferentiate between models. Both optical (e.g. ER 0047-2808) and 
radio (e.g. MG075 1-1-27 16) lenses can fall into this category. 

An alternate strategy for defining in the case of a "thin an- 
nulus" image is to define a mask in the image plane and only calcu- 
late over pixels in the mask. The mask should cover the image 
regions (plus a little extra so that incorrect models can be appro- 
priately penalised) and should ensure that all images of a source 
pixel are included in the mask. In this case, the mask will cover 
pixels both inside and outside the tangential caustic, so Equation 
cannot be used. Instead we define 

V = A^niask — A^sourcc ^ m (8) 

where A^mask is the number of beam areas inside the mask (or pix- 
els for CCD data), A^'sourcc is the number of pixels in the source 
plane which map to the mask region and m is the number of lens 
model parameters. 

There are two issues that must be raised relating to the number 
of degrees of freedom used in the source plane. 

Firstly: the entropy constraint introduces correlations in the 
source plane. How does this affect the number of degrees of free- 
dom in the source? In the context of the lens modelling problem, 
this is not a serious issue. The entropy constraint works on the over- 
all value of J, so pixels are connected only loosely via their con- 
tribution to J. The gradient of the entropy, VSki, is applied to in- 
dividual pixels so there is no correlation between pixels introduced 
by this operation. When the algorithm is used to find the best pos- 
sible x^, the source model is determined entirely by the data with 
the exception that source pixels cannot become <^ A (or negative). 

Secondly: do source pixels with small values (e.g. less than 
the noise) contain degrees of freedom? It is possible to reconstruct 
sources having a significant number of pixels which are effectively 
zero, in the sense that they are well below the noise level in the im- 
age. These source pixels contribute (by definition) almost nothing 
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to the . Yet in our definition {Sj for the number of degrees of free- 
dom, these pixels are counted as free parameters which artificially 
lowers the target . Instead, a threshold, T, can be defined, be- 
low which all pixels are set to zero. Two models are different (with 
99% confidence) if their values differ by 6.63, so T should be 
defined such that the overall change to is ^ 6.63 after setting 
source pixels to zero. In this way only, pixels which contribute sig- 
nificantly to the model image (therefore x^) are counted as degrees 
of freedom. This idea is demonstrated in f|3| 

For radio data, where the noise between neighbouring pixels 
in the image plane is correlated, problems defining a target (cal- 
culated in the image plane) remain. The use of a pixelised source 
plane and the entropy constraint makes Equation an impossible 
target to reach, even for the correct lens model. Our tests found this 
to be true even with the entropy constraint disabled. The problems 
arise because the basic iterative technique (using gradients of x^ 
and the mapping matrix) is fundamentally different to a CLEAN- 
based approach where a resolved source is constructed from many 
point sources. The entropy constraint means that data values that 
are negative also contribute more to the x^ than they would using 
CLEAN components. Since regions of negative brightness close to 
the real sources are a common artifact of CLEANed radio maps, 
this effect must be compensated for somehow. 

Lensview does not yet work in visibility space, so the prob- 
lem of defining a targ et x^ for radio data remains u nsolved for now 
although the work of iKemball & Martinsekl i2005h in determining 
the uncertainties of pixels in radio maps shows promise for the fu- 
ture. This means that, for now, it is not possible to say whether a 
particular lens model is an acceptable fit to the data using radio im- 
ages. It is shown in ^ however, that lens models can be compared 
reliably using LENSVIEW with radio images and that lens model 
parameters are reliably recovered by the software. For optical data, 
where errors in neighbouring pixels are uncorrelated, these prob- 
lems do not exist. Using optical data, lens models can be quantita- 
tively evaluated by their ability to model the data. 

2.4 Reverse Projection 

WNK94's equation 8 showed how to calculate dJ/dSi (VCti and 
VSfci in our notation) in one dimension. This equation incorpo- 
rates the beam and mapping matrix for reverse-projecting VCij as 
expected. However, it is not completely correct for the reason dis- 
cussed below. 

The beam must be modified before it is applied as part of 
VCij in the image plane. We consider for a moment only the im- 



age plane with beam-convolved model image Mi, 



Br, 



denotes convolution. Before 
reverse-projecting into the source plane, we must calculate VCij = 
dC/dlij. Equation 8 of SB84 showed (using one-dimensional no- 
tation with no loss of generality) that the transpose of the matrix 
representing the beam must be used in calculating VCi-,, which 
was not used in WKN94. It can be shown, after some algebra, that 
the equivalent transformation for a two-dimensional beam operat- 
ing on two-dimensional images is to rotate the beam image by 180 
degrees. Thus, WNK94 Equation 8 is correct for beams that are 
symmetric under a 180° rotation (e.g. in radio data), but will not be 
correct for asymmetric beams. 
We can write 
dC 
dfki 



dC dlij 



dhj dfki 



(9) 
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Figure 1. Example tiling of the image plane onto a section of the source 
plane. The centre of the source plane is at x=0,y=0 (not shown). Projected 
pixels from two distinct regions of the image plane are shown (shown in 
parentheses, or by row/column number). Pixels around (49,49) come from 
very close to the centre of the lens, whereas pixels around (83,87) are at the 
edge of the multiply imaged region of the lens. 



dC 

dlii 



,(A: 



Mi, 



(10) 



is the gradient of x including convolution with the rotated beam, 
B^, and 



Wii 



ijkl 



follows from Q. So, to reverse-project VCij we calculate 

^CijWijki 



(11) 



(12) 



where 



where iTOT.ij = Sfc'!' ^ijk'V = ^imagc/Mij IS the total area 
mapped to by an image pixel. 

An understanding of the subtleties of reverse projection can be 
gained using an example. Consider a simple elliptical isothermal- 
style lens. A model is made using a 100 x 100 image plane and 
a 60 X 60 source plane. The centre of the lens in the image plane 
is pixel (52, 52)i. (Image plane coordinates are denoted with sub- 
script i and source plane with subscript s to avoid confusion). In 
Fig- a section of the source plane is shown with pixels from the 
image plane projected onto it. There are two distinct regions of 
the image plane that map to the region in the source plane around 
(19,19)3. The block of projected pixels from around (49,49)i 
originates from very close to the lens centre and are highly de- 
magnified. The block of pixels from around (83, 87)i originates 
from just inside the boundary of the multiply imaged region and 
are slightly magnified. Source pixel (19, 19)s contributes to (and 
hence will take contributions from in reverse-projecting VCij) pix- 
els (50,50),, (49,49),, (49,50),, (82,86),, (83,86),, (84,86),, 
(82, 87)i and (83,87)i. One can see that pixels (83,87)i and 
(50, 50)i are the two pixels with roughly the largest weight and 
VCfcj would take a roughly equal contribution from both (from 
Wijki) with lesser contributions from the remainder. However, 
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pixel (50, 50)i is strongly demagnified and source pixel (19, 19)s 
has only a small fraction of the total area in the source plane to 
which it maps. The correct scheme then, is to weight each image 
pixel's contribution by the fraction of the total projected area that 
a source pixel contains. This is the same as weighting an image 
pixel's contribution by its magnification. This factor explicitly ap- 
pears in the denominator of J12t . What about source pixels in re- 
gions of high magnification? In such a region, there are several (2 
or more) distinct parts of the image plane mapping to the same part 
of the source. In addition, several neighbouring pixels within each 
image region are mapped inside the same source pixel or small re- 
gion. Again, there is a problem because there are many image pixels 
all contributing to the source, none of which are down-weighted. 
In this case, the total contribution to the source pixel will be too 
large, by a factor that is equal to the magnification at that location 
in the source plane. To solve this problem, the overall contribution 
to a source pixel must be down- weighted by ^ki- This factor is ac- 
counted for by the definition of Wtjki- 

We note also that in testing the reverse projection process we 
tried several alternate definitions of Equation <12> . We found that 
some work well except in regions close to the boundary of the mul- 
tiply imaged region, such as Fig.Q In that case we found that flux 
was suppressed in all source pixels to which (50, 50)i mapped (or 
the equivalent region). The algorithm tried to compensate for this 
by raising the value of surrounding pixels, causing sudden (unre- 
alistic) jumps in the pixel values which were reminiscent of the 
"glitches" reported in WNK94. It is possible that the glitches were 
caused by an incorrect reverse-projection process. We find no such 
problems using reverse-projection defined by <I2t . 

2.5 The Skilling & Bryan MEM algorithm 

SB 84 provide several key insights on using iterative methods with 
MEM, which are summarised here. 

The bulk of the computational cost in calculating Afki comes 
from convolution (via FFTs) with the beam to create the model im- 
age and convolution with the rotated beam to create Vdj (plus 
forward and reverse projection for lens modelling). If the conju- 
gate gradient method is used for each inner loop iteration, a sin- 
gle search vector A/^; is calculated then a one-dimensional line 
minimisation is performed to find the scale factor (e) such that 
/ -I- eAfki minimises J for that iteration. SB84 noted that by 
combining VCki and VSti into a single vector, useful informa- 
tion about the properties of J is being discarded. In addition, be- 
cause A is used to combine VCki and VSki before the line min- 
imisation is performed, it cannot be adjusted as part of the process. 
SB 84 found a more general and flexible model can be generated, at 
negligible extra computational cost, by using a subspace of search 
vectors based on VCki and VSki- 

Lensview replaces the conjugate gradient based inner loop 
of LensMEM with SB84's method. The basic idea is summarised 
here and the interested reader is referred to SB84 for details. For 
each iteration of the inner loop, a subspace is constructed spanned 
by three search vectors: 
ei =/(VSfcO, 
62 =f{\/Cki) and 

ea =f{VVCki)f{VSki)/\VSki\-f{VVCki)f(VCki)/\VCki\. 
Note that both of the original search directions, (VSki, VCki ), are 
weighted by /, the "entropy metric", (pixel by pixel multiplication, 
not a dot product) which is necessary and useful, but will not be dis- 
cussed here. VVCfc; is simply the reverse-projected per-pixel vari- 
ance image after convolution with the beam and the transformed 



beam and is constant throughout the inner loop. The subspace is 
simply a model of the local properties of J and depends only on 
and A. 

The goal of the process now is to find x'' such that /+ 
minimises J in the subspace, subject to some sensible constraints 
on how much the source can be changed in a single iteration. In 
addition, because VSti and VCki are not combined into a single 
search direction, A can be adjusted as part of the same process. 
Indeed, for each iteration of the Skilling & Bryan algorithm, the 
minimum in J (within the subspace) always represents the maxi- 
mum entropy solution for that iteration. In principle, and A can 
be solved for algebraically, however constraints on how much the 
source can be sensibly changed in an iteration mean a more compli- 
cated process is required. The actual process for determining A is 
via the "a-chop" algorithm described in SB84, which does not use 
any projections or convolutions. Once A is determined (per itera- 
tion), x'' is calculated, a new source generated and the next iteration 
begins. 

The definition of entropy has also been adopted from SB 84 
(Equation |5j, which incorporates a default source value, A. This 
is the value that a source pixel will take if there are no data con- 
straining it. Introducing A makes S a pseudo-dimensionless quan- 
tity in the sense that any source pixel which has a value equal to 
A will contribute zero to the entropy, regardless of the units of 
measurement of fki. In contrast, the definition used by WKN96, 
S — fki Infki, does depend on the units of measurement 
whereas the statistic does not. In practice this unit difference 
can be absorbed into the value of A, however pixels with values 
<^ 1 are strongly prohibited by the WKN96 definition of entropy 
where dS/dfki ~ —1 — In fki - The SB84 definition of entropy has 
dS/dfki = In j4 — In fki which tolerates values that are not <^ A. 

The need to define A, the default source value, highlights the 
importance of understanding the noise characteristics of the data. 
This point has been stressed many times in the context of the MEM 
and is just as important for LENSVIEW. The noise properties of the 
data affect the results in the calculation of C, the default source 
value and in the calculation of VCij and WCki- A suitable value 
for A is roughly 1/2 the value of the background noise in the image. 
The results of modelling are not very sensitive to A as long as it is 
the same magnitude as the noise. 

Lensview is implemented such that it can use the Lens- 
MEM-style conjugate gradient approach or the SB84 algorithm. 
Extensive testing has shown that the SB 84 algorithm finds better 
sources* (both smoother and better x^) and is much faster than the 
conjugate gradient approach. This is because the multiple forward 
projections to find e in each iteration of LensMEM's inner loop are 
removed. The conjugate gradient approach is useful if one wants to 
disable the entropy constraint by setting A = 0. 



3 TEST CASES 

There are two main goals when testing the software. Firstly: can 
the software recover the lens model parameters of a known lens, 
and with what accuracy? Secondly, is the source reconstruction ac- 
curate? In this section, the software is applied to simulated optical 
and radio data with these two goals in mind. While working to- 
wards the high-level goal of answering the two questions above, we 

* This is due mostly to the difference in definition of entropy and to a lesser 
extent the definition of search directions (V/ vs e^). If LENSVIEW used 
WKN96's definition of entropy, the results should be the same in principle. 



© 2002 RAS, MNRAS 0Q0.[T1|22| 



6 R. B. Wayth and R. L. Webster 



must also consider the issues of source plane pixel size, the overall 
size of the source plane and the number of degrees of freedom used 
in the source model. 

In both tests, the data is processed in image space. For optical 
data, the analysis is straightforward using the data and a pixel vari- 
ance image. For the radio data, the image is not directly measured 
by the telescope, hence the analysis must begin with understanding 
what effects, if any, are introduced by using processed images. 

The focus of each section will be different because the typical 
lensing scenario is quite different between (say) a lensed high red- 
shift galaxy (typically ~ 0.1 arcsec in size) and a lensed radio lobe 
(~ 1 arcsec in size). 

For a lensed galaxy, the object is of comparable size to the 
tangential caustic and usually highly magnified into a thin ring or 
series of arcs (e.g. the lenses ER0047-2808 and the near-IR images 
of QSO host galaxies). The high magnification reveals structure 
in the source that would be unresolved if not for the lensing and 
is manifested through brightness variations in the ring/arcs. This 
structure must be accounted for in the source model so the spatial 
resolution of the source plane must be higher than the measured 
image. The overall size of the source plane, however, is small. 

For a lensed radio lobe, the source is much larger than the 
tangential caustic, hence the overall magnification is modest. The 
object is lensed into a ring with surface brightness that varies, but 
is relatively constant around the tangential critical line. Hence it is 
not likely that source pixels need to have higher spatial resolution 
than image pixels (which are themselves typically 1/3 of the beam 
width). In addition, much of the image will be found away from 
the high magnification region but still multiply imaged. The source 
plane in this case needs to be large. 



3.1 A test with simulated optical data 

In this section, Lensview is tested with simulated optical images. 
The images are designed to test the discriminating power of the 
data for images with significantly different sizes. The issues that 
specifically need to be addressed are: the optimum size of pixels in 
the source plane, the accuracy of lens model parameters, the dis- 
tinguishing power of the image between different (non-isothermal) 
models and the overall goodness-of-fit of any particular model. 

Three model images were created using a two-component 
source. The source and images are shown in Fig. |2| The source 
is designed to have two regions of substantially different size and 
different peak brightness. The source was offset from the lens 
centre and projected into an image with 0.05 arcsec pixels us- 
ing a softened power-law elliptical mass distribution (SPEMD) 
( Kassiola & Kovner 1993; Barkana 1998) lens model with non- 
isothermal mass slope. Each of the model images has the same 
lensing properties, it is only the relative size of the source that has 
changed. The SPEMD is defined by surface mass density 



Ratio X 



Comment 



/ 2 , 2 I I 2n~/3/2 

K = K,o{x q + y /q + r^) 



(13) 



where kq is the central density, q is the axis ratio, Vc is the central 
core radius, f3 is the mass slope parameter and the major axis of the 
system is rotated by 9 degrees. The mass slope allows the SPEMD 
to model a range of mass profiles from constant density (/3 = 0) to 
isothermal (/3 — 1) to a modified Hubble profile (/3 — 2). In prac- 
tise, we scaled kq such that b — kq^^ (^^) ^ ^^^^ ^ single 
parameter, b, is used for the mass scale factor which represents the 
Einstein radius of the lens. The parameters used to create the sim- 



1 1266.6 Source components not resolved. 

Poor image reconstruction. 

2 1 146.5 Barely resolved source. 

Improved image reconstruction. 

3 1086.9 Resolved bright regions in source. 

4 1055.8 Source becoming noisy 

5 1024.6 Source very noisy 

Table 1. Best that could be reached using different image-to-source pixel 
size ratios for simulations of optical image 3. The was evaluated over a 
mask of 1234 pixels surrounding the image. 



ulated image were b = 1.2 arcsec, q — 0.82, Tc — 0.05 arcsec, 
/3 = 1.16 and 6 = 68°. 

A PSF was generated for the HST ACS instrument using the 
F555W filter with the TinyTim <Krisllll993h software. This PSF 
was convolved with the simulated image. Next, random noise was 
added to the image including contributions from CCD read noise, 
sky and signal shot noise. (The effects of lensing galaxy subtraction 
have been ignored here and it is assumed that the lensing galaxy has 
been subtracted perfectly. In addition, it is assumed that the lens 
galaxy light does not add significant noise to pixels in the region of 
the image.) The overall signal-to-noise ratio was chosen to be sim- 
ilar to ER0047-2808 with the peak signal in the simulated image 
approximately 10 times the background noise level. The sources, 
PSF-convolved images and noisy images are shown in Fig.|2| The 
model images are labelled 1-3 as shown. 



3.1.1 The optimal source pixel size 

When reconstructing a source, pixels in the source plane need only 
be small enough to explain the smallest features in the image. In 
this section we investigate the best pixel size for image 3 because it 
has the finest detail. 

An estimate of the smallest required pixel resolution can be 
made by examining the image data. The arc crossing the tangential 
critical line has two bright regions on each end. These regions are 
resolved into small arcs about 3 pixels long. Given that these im- 
ages are in a high (~ 10) magnification region and all the magni- 
fication in isothermal-style lenses comes from tangential stretching 
of the image, the size of the structure in the source must be ~ 3/10 
the size of an image pixel. Thus, a first guess for a suitable source 
pixel size is 0.3 times the image pixel size, i.e. an image-to-source 
pixel resolution ratio of 3.33. 

To test both the ability of the software to reproduce the data 
and the quality of the source reconstruction, image 3 was analysed 
with source plane pixels I, 1/2, 1/3, 1/4 and 1/5 the size of the 
image pixels (i.e image-to-source pixel ratio of 1, 2, 3, 4 and 5 
respectively). A mask of 1234 pixels surrounding the arcs was de- 
fined and the image was processed using the known (correct) lens 
model parameters. The source plane overall size was chosen to fit 
the reconstructed source with a little extra room. 

Results of the modelling are shown in TableQand graphically 
in Fig. |3| The quality of the model image improves dramatically 
for a pixel ratio of 2 (compared to I) and again for a pixel ratio of 
3. For pixel ratios of 3 and above there are no discernible residuals 
between the data and model image. The meaning of the values 
depends on the number of degrees of freedom in the source — a 
topic that is discussed later. However, it is clear that a pixel ratio 
of 1 produces a poor model image. For higher resolution source 
planes, the source reconstruction becomes very noisy because dis- 
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Figure 2. The source and model images (noiseless and noisy) used in the optical data simulation. The source is shown with the caustic to highlight its position 
relative to the lens centre. The noise-free image is shown with the tangential critical line. All images have a linear greyscale stretch. 
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Figure 3. Results for reconstructions of simulated optical data (test image 3) with different source pixel sizes. The source planes have been scaled to the same 
angular size. All images have a linear greyscale stretch. The image on the left in the row "difference images" is the actual image, including noise, used in the 
modelling. 
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Figure 4. Reconstructed sources for the simulated optical data with differ- 
ently sized source planes. The sources are 10x10, 15x15, 20x20, 25x25 and 
30x30 from left to right. For all images, the image-to-source pixel size ratio 
is 3. 



Parameter 


True value 


mean 


std dev 


Skewness 


b (arcsec) 


1.2 


1.200 


0.0017 


0.08 


q 


0.82 


0.80 


0.02 


1.73 


9 (degrees) 


68 


68.03 


0.42 


-0.19 


/3 


1.16 


1.22 


0.02 


1.61 






1142 


45 


0.04 



Table 3. Statistics of the fitted parameters for 400 realisations of the simu- 
lated optical source. 



tinct regions of the source are able to model individual noisy pixels 
in the image. Hence they are affected by the noise in the image. At 
medium resolution (ratio of 2-3), the source is detailed enough to 
be able to reconstruct the image but the pixel mapping overlap from 
neighbouring pixels helps to smooth out the source. 

From this qualitative analysis, it is clear that an image-to- 
source pixel scale ratio of at least 2 is required to reproduce the 
data and a ratio of 3 is optimal for this data. This agrees well with 
the initial guess for a ratio of 3.33. Also shown in Fig.|3|are the 
source reconstructions for the noiseless (but still PSF-convolved) 
image. In that case, the reconstructions match the original source 
very well and are limited only by the finite resolution of the image. 

3.1.2 Model goodness-of-fit and source degrees of freedom 

Using Equation (Sj, an acceptable range for can be defined for 
the simulated data. The equation relates the degrees of freedom 
used in the source plane directly to the target . In this section, 
this idea is explored further. 

Using an image-to-source pixel size ratio of 3 and a variety 
of source plane sizes, the best possible the model could achieve 
(subject to the entropy constraint, of course) was found, again for 
test image 3. In each case the source plane is positioned so the re- 
constructed source is as close to the centre as possible to minimise 
any potential edge effects. The resulting reconstructed sources are 
shown in Fig.|4| 

It is clear from these images that once the source plane is large 
enough (approximately 15 x 15 pixels in this case), most of the 
remaining pixels in the image do not contribute anything to the 
model. Using the pixel thresholding scheme presented in i|2.3l the 
simulated data was analysed again with pixel ratio 3 and the 1234 
pixel mask using a 99% confidence threshold for T. Hence, was 
permitted to change by up to 6.6. Table|2|shows the best that can 
be reached using different sized source planes and the interpretation 
of the model goodness-of-fit both with and without the thresholding 
scheme. 

The table shows that large source planes have a minimal ef- 
fect on improving the model but they do help to clearly identify the 
location and extent of the source. (For very large source planes, the 
pixels will map outside the pixel mask so by definition do not con- 
tribute to .) Hence, the simple definition {5} for source degrees 
of freedom is oversimplified for all source planes which have many 
low-value pixels. The threshold-based definition means the results 
do not depend on the particular size of the source plane (unless it is 
too small), which is the result one intuitively expects. 

Using the threshold based definition for v, 

V = Nma.sk - (iVsourcOT) - m (14) 

it is clear from Table |2| that all reasonably sized source planes are 
able to produce a model which is statistically acceptable. Only the 
10 X 10 model was not within 99% confidence because it was too 



small to contain the source. For large source planes, regions of the 
plane that map outside the image mask are not affected by the data 
hence will become the default pixel value, A. These pixels then do 
not affect the fit at all because they have zero contribution to J and 

Given that the noise properties of the image and lens model are 
known exactly, it is very reassuring that the thresholding scheme 
produces the expected results. 

3.1.3 Errors on fitted parameters 

The focus of testing thus far has been on the properties of the source 
and source plane using a fixed (correct) lens model. Obviously the 
lens model parameters are not known with real data and the major 
outcome of modelling is a set of measured lens model parameters. 
Thus, the next task is to quantify the uncertainty in the lens model 
parameters themselves. 

In a linear modelling experiment, the value of x^ should (lo- 
cally) form a quadratic hyper-surface around the minimum. In such 
a case, the uncertainties on fitted parameters can be estimated 
dire ctly from the properties of the surface [e.g. section 15.6 of 
|PressetalJjl992h l. 

An image generated through gravitational lensing is highly 
non-linear in the lens model parameters. Hence another technique 
such as Monte-Carlo simulations must be used to estimate the 
uncertainties on fitted parameters. For very small changes in the 
model parameters, however, the result might be 'locally linear' so 
it is worthwhile examining the properties of the x^ surface anyway. 
If the x^ surface locally approximates a quadratic, then quantify- 
ing a parameter's uncertainty is a much easier task than having to 
do extensive simulations. In this section, a Monte-Carlo simulation 
is performed on the simulated optical lens to measure the true un- 
certainty in fitted parameters. These results are then compared to 
the values derived from the x^ surface. 

Using the same simulated optical image and noise as before, 
400 datasets of the simulated optical source were generated. For 
each simulation the best fittin g parameters wer e found with the 
downhill simplex method (e.g. lPress et alJ^'992^ using a SPEMD 
model with parameters (b, q, 9, /3) free to vary. The lens centre and 
core radius were fixed because the lens centre is usually well deter- 
mined for an optical lens and there are no central images, so there is 
no point investigating a core. A 15 x 15 source plane was used with 
default source pixel value (the A parameter from ij2j fixed equal to 
the background noise level in the image. 

The mean, standard deviation and skewness of the best-fit sim- 
ulation parameters are shown in Table |3| The table shows that all 
of the lens model parameters have been well recovered although in 
the case of the axis ratio and mass slope, the best-fitting parameters 
are not distributed normally. The tail in the distribution increases 
the standard deviation (hence uncertainty) of a measurement. (A 
skewness of 1 makes the distribution approximately twice as wide 
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1.0 
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Table 2. Results for differently sized source planes that have source pixels < T set to zero. Column 1 gives the size of the source plane in pixels. Column 2 is 
the best attainable with that source plane. Column 3 is the target using equation |8). Column 4 lists whether the is within a 99% acceptability range. 
Column 5 lists the threshold value used for the source plane. Column 6 lists the number of pixels in the source that are above this threshold. Column 7 lists the 
revised target using equation I14i . Column 8 lists the for the source with pixels below T set to zero. Column 9 lists the acceptability of the new model 
with the revised target. 
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Figure 5. Best-fit parameters over all the realisations of the simulated optical source. 



compared to a normal distribution.) Interpreting the standard devia- 
tion value as 'Icr' error, the mean of the best-fitting parameters are 
all within lu of the true value, with the exception of 13. We note 
that of all the model parameters, only /3 seemed to be affected by 
the starting position of the simplex search and even then the results 
only varied by /?± ~0.02. Noteworthy is the correlation between 
(3 and q as revealed in the Fig. |5| This correlation is related to the 
ellipticity/m ass-power-law degeneracy found in position-only QSO 
lens models iWambseanss & Paczvnskll994) . 

Using a single simulation, parameter space was searched 
around the minimum to make a map of the surface. Shown in 
Fig |6(a)| are marginalised contours of for various combinations 
of parameters. The contours are 2.3, 4.6, 9.2 and 18.4 above the 
minimum which correspond to the 68%, 90%, 99% and 99.99% 
confidence intervals for the two parameters shown. 

Examining the first two plots in Fig |6(a)| it is clear that the 
mass scale, 6, is well constrained and that the 68% confidence re- 
gion (approx 1.200 to 1.204) corresponds well to an uncertainty of 
±0.002 suggested by the Monte-Carlo simulations. The third fig- 
ure of Fig |6(a)| shows a large uncertainty in q and /3- far more than 
is suggested by the Monte-Carlo simulations. In this case the 13 — q 
degeneracy forms a valley through space, although it is reassur- 
ing that the simplex method used in the Monte-Carlo simulations 
has found its way down in most cases. Again, if we interpret the 
la errors from the 68% contour and from the Monte-Carlo simu- 
lations to be /3 ± 0.05 then, the true value is less than 2a from the 
minimum at /3 = 1.23, q — 0.79. The implication for measuring 
errors is that the long valley of 99.99% uncertainty does not seem 
to affect the results. We conclude then, that the surface is a good 
measure for the statistical uncertainty of the lens model parameters. 

As an illustration of the accuracy to which b can be measured 
if f3 is fixed, a slice (not marginalised) through the surface is 
shown in Fig. |6(b)| (LHS). In this case, the 68% uncertainty is 0.004 



arcsec. Similarly, a slice through (3 — q space (RHS) shows a well 
localised minimum. Hence, if extra constraints can be added to a 
lens model through (say) morphological or velocity dispersion in- 
formation, the resulting lens model parameters will be more tightly 
constrained. 



3.1.4 Quantitative comparisons between lens models 

Because the properties of lensed images are determined mostly by 
the lens mass profile in the region of the images and the 'correct' 
lens model has been used in these simulations, it is not surprising 
that the mass slope has been recovered well. This leads to the ques- 
tion: are there other lens models with a similar mass slope in the 
region of the images, but quite different elsewhere, that can fit the 
data? This question is the topic of this section. 

While there are virtually limitless combinations of compo- 
nents that could comprise a mass model, the main models of in- 
terest here are the singular isothermal ellipsoid in external shear 
(SIE-1-7) and the constant M/L model. These models are motivated 
by the fact that the image positions in almost all QSO lenses can be 
fitted with a SIE-1-7 model ( Keeton et al. 1997) and that a power-law 
mass distribution (SPEMD) has the potential to mimic the effects 
of a constant M/L profile under some circumstances (if the mass 
profile slope is locally constant at the image radius). The SIE is the 
same as a SPEMD with (3 — 1. The external shear is parametrised 
by its magnitude, 7, and position angle, 6.^ . 

Testing the SIE-1-7 model is important because external shears 
can be used to improve lens models which otherwise provide poor 
fits to the data. If an external shear can mimic the effect of a non- 
isothermal mass distribution, then there is little hope that accurate 
mass profile slopes (i.e. f3) can be measured using image data alone. 

The case of the constant M/L model is much harder to test 
because there is no photometric data associated with the simulated 
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Figure 6. Contours of for different combinations of model parameters. Contours increase from the minimum by 2.3, 4.6, 9.2 and 18.4. Crosses show the 
parameter values used to create the data. 



optical source. Typically, the photometric properties of the lensing 
galaxy are quite well determined, hence a test for the acceptability 
of a constant M/L model is fairly straightforward: just build the 
lens model based on the photometric parameters and vary the M/L 
to find the best fit for an absolute goodness-of-fit measurement. 
(Whether the best fit M/L is realistic is another matter). For this test 
case, it will have to suffice to test whether a constant M/L model 
can mimic a power-law model rather than the other way around. 

The con stant M/L mass distribution is modelled with a Sersic 
jSersicll968h profile. The Sersic profile is defined by 

K(r) = Koexp-^<")(''/''=''^" (15) 

where is the central surface density in units of the lensing crit- 
ical density (Ecrit ~ (? Ds / A-KGDdDds), fa is the scale radius 
and n is the shape parameter. The profile can also be elliptical with 
axis ratio q, and position angle 8. The Sersic profile incorporates 
the de Vaucouleurs (n = 4), exponential (n = 1) and Gaussian 
(n — 1/2) models with higher n producing more cuspy central 
peaks. In this test, n is fixed to 3 based on the photo metric pro per- 
ties of the lensing galaxy in ER 0047-280 8 ( Wavth e t all2005h and 
for elliptical galaxies in general jBalcells e t al. 2003). 

A constant M/L model that has the best chance of mimicking 
a power-law must have a mass slope which is locally constant and 
equal to /3. Shown in Fig.0are plots of deflection angle vs radius 



for circularly symmetric Sersic profiles with n — 3. Also shown is 
the deflection angle vs. radius of a power-law model with /3 = 1.16. 
This figure shows that for n = 3, Sersic scale lengths, r^, between 
approximately 2 and 3 produce deflections that are very similar to 
the power-law model in the region of interest. Hence, our tests for 
the constant M/L model were limited to n = 3 and a range of 
between 2 and 3. 

Using the same simulated data as the previous section, param- 
eter space was searched for SIE, SIE+7 and constant M/L models. 
The results of the search are shown in Table|4] 

The first result from this test is that adding an external shear 
to an SIE model cannot improve the fit compared to a plain SIE 
model. Hence, adding an external shear cannot help improve the 
model for this type of lensed image. This result is consistent with 
the fact that an external shear distorts the critical line away from 
being an ellipse. Since the image has flux at various points near the 
critical line, the image is sensitive to such distortions and can rule 
them out. 

The second result from the test is that the constant M/L model 
can produce a result which is better than the SIE. All of these mod- 
els produce statistically acceptable results although the difference 
in reveals that the models are different with > 99% confidence. 
Hence, if this was real data the SPEMD model would be deemed 
the best model. 
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Model Best Parameters 

SIE 1138 6 = 1.225,9 = 0.854,0 = 67.6 

SIE + 7 1138* 6= 1.225, (j = 0.854, = 67.6, 7 = 

Constant M/L 1121 kq = 57.9, 5 = 0.85, 6» = 67 

SPEMD 1107 b = 1.200, g = 0.809, e = 67.7, /3 = 1.16 (fixed) 

Table 4. Best-fitting models for the simulated optical data with incorrect lens models. *Note: No models with improved could be found using an SIE with 
external shear. 



2.5 



2.0 -> 



1 .5 



1.0 



0.5 



0.0 L 
0.0 



power law M = r ^-^^ 
Sersic. n-3 r,-l 
Sersic. n-3 i\-2 
Sersic. n=3 T^=4 



1.0 1.5 

Radius (arcsec) 



2.0 



Figure 7. Relative deflection magnitude for Sersic profiles with different 
shape parameters in comparison to a power law model. 



These tests have shown, therefore, that this type of optical im- 
age can indeed distinguish between lens models although the differ- 
ences (in terms of x^) are quite small. Statistically speaking, both 
the SIE and constant M/L models are acceptable in terms of the 
data alone. 



3.2 A test with simulated radio interferometer data 

In this section, a series of tests is performed with simulated radio 
data. To make the tests realistic, the lobe must be 'observed' with 
a simulated telescope then an image generated for the software to 
analyse. Typically, the image is generated with the CLEAN algo- 
rithm tHoebom 1974;.qark^,1 98Q.:.Schwab,1984: Steer et al..l984i) . 
see also jCornwell e t aljl999l) . Differences between the real image 
and the CLEANed model can arise because of incomplete (u, v) 
coverage and/or gain and phase errors introduced in the observing 
process. To distinguish between these effects, it is necessary to cre- 
ate simulated data for different length observations both with and 
without errors. 

Simulated radio data were created by taking the unlensed ra- 
dio lobe from MG1549-I-3047 (shown on the left in Fig.[n} as the 
'true' sky brightness (as seen by a noiseless telescope with infi- 
nite resolution). The lobe was corrected for any negative values 
then projected through a PIEP lens (defined by the lensing potential 
ij) = &(r^ + (l-e)a::^ + (l + e)y^)^''^ with Einstein radius, b, ellip- 
ticity, e, orientation angle, 9 and core-radius, Vc) with b — 1.35", 
e = 0.065, 9 = 0° and rc = 0.^ The location of the lobe (relative 



For rc 0, b is not the Einstein radius, it is simply a mass scale factor 



r 




Arcsec West 



Figure 8. The simulated lensed lobe 'true' image. Contours increase by \/2 
from 0. 1 mJy beam" ^ . 



to the lens) was chosen such that the bright hot-spot was close to 
the inside edge of the radial caustic while the fainter, diffuse tail 
of the lobe passed near and through the centre of the lens. It was 
assigned a position in the sky near MG1549-H3047. The simulated 
lensed lobe is shown in Fig.|8| 

Using the MiRIAD tasks 'uvgen' and 'uvmodel', visibilities 
were generated for observations of the simulated lobe with the 
VLA A array at 8.4 GHz. Observations were generated both with 
and without realistic gain and phase noise for 1-hour and 3-hour 
integrations. Observation time was spread approximately evenly 
around HA = 0. Plots of the {u, v) coverage for both integration 
times are shown in Fig.|9| 



3.2.1 Image reconstruction with ideal and realistic data 

Images from interferometers with incomplete (m, v) coverage must 
be reconstructed with some algorit hm like CLEAN. This introduces 
known artifacts into the image Isee lComwell et alj |T999) for a de- 
tailed discussion]. In the context of lens modelling, it is the effects 
of these artifacts on the lens model that are of interest. 

The images were processed with the MiRIAD tasks 'in- 
vert' 'clean', 'selfcal' and 'restor' using 10000 iterations, natural 
weighting, only positive clean components and a small (0.02) loop 
gain. The field size was 512 x 512 pixels of size 0.05 arcsec and the 
beam FWHM 0.28 x 0.27 arcsec. For the simulated observations 
that contained gain/phase errors, the 'selfcal' task was used to cor- 
rect the complex gains using self-calibration. Calibration solutions 
for the VLA are over-constrained, therefore well determined. In all 
cases of self calibration, the model converged in a single iteration. 



which is larger than the actual critical radius of the lens. Nevertheless, it 
will still be referred to as the 'Einstein radius' or the 'critical radius'. 
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(b) 3 liour integration 

Figure 9. {u, v) coverage for simulated 1 hour (top) and 3 hour (bottom) 
observations of the lensed lobe. 

Shown in Fies. llOl and llll are the CLEANed images for the sim- 
ulated observations along with the 'true' image convolved with a 
Gaussian beam of the same size as the CLEAN restoring beam (i.e. 
the ideal reconstruction). 

On the contour plots, the differences between the images are 
most obvious in the region of fainter, diffuse emission north and 
west of the inner image. However, it is not the faint, low magni- 
fication emission that provides the constraints for the lens model, 
it is the bright arcs and ring. In Fig. 1121 the differences [i.e. the 
ideal image of Fig. 1101 (top) minus the reconstructed image] be- 
tween the beam-convolved ('true') image and the CLEANed im- 
ages are shown as images, not contours. The artifacts generated by 
CLEAN are more apparent in these images. The large arc has a 
thin shadow around most of the bright region and the inner bright 
image is misshapen. Because the maximum/minimum of the differ- 
ence between the bright regions in the images is on the same side 
for both images, the effect is that the image appears to be shifted. 
However, this effect is a coincidence. There has been no transla- 
tion of the CLEANed image. The effect is likely a consequence 
of the finite pixelised grids generated by Fourier transforming and 
CLEANing the data. CLEAN components are chosen from a grid, 
so there must always be a small (fraction of a pixel) difference in 
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Figure 10. Top: The simulated lensed lobe convolved with a Gaussian 
restoring beam which is the same size as that generated from the 1 hour 
simulated observation. This is what a perfectly CLEANed image should 
look like. Bottom: The CLEANed image from a 3 hour simulated observa- 
tion of the lensed lobe with no observation errors. Contours increase by \/2 
from 0. 1 mJy beam^ ^ . 



position between the peak flux suggested by the data and the loca- 
tion of the corresponding CLEAN component. It is apparent from 
these images that the main source of error in the image reconstruc- 
tion comes from incomplete {u,v) coverage and CLEANing, not 
observational/instrumental errors. 

Given that most of the information about the lens model will 
come from the bright regions, the outcome of the fitting process is 
likely to be affected. In particular the apparent translation of the 
reconstructed images will certainly affect the determination of the 
lens centre. It is also worth noting that the self-calibration process 
can lose absolute astrometric accuracy in radio data so radio im- 
ages should be aligned with optical whenever possible to avoid this 
problem. 



3.2.2 The effects of CLEANing on lens modelling 

To quantify the effects on the lens model of the artifacts discussed 
in the previous section, a series of Monte-Carlo simulations was 
performed. One hundred simulated observations of the simulated 
lensed lobe from the VLA A array at 8.4GHz were generated with 
realistic system noise, gain/phase errors and an integration time of 
1 hour. The observations were processed by an automated script 
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Figure 11. CLEAN reconstructions of the simulated observations of tlie 
lensed lobe. Top: The image from a 1 hour simulated observation with 
no errors. Bottom: The image from a 1 hour simulated observation with 
gain/phase eiTors. Contours increase by \/2 from 0.1 mjy beam~^. 




Figure 12. Differences between the beam-convolved true lens image and 
the CLEAN reconstructions. Left: the 3 hour observation with no errors. 
Centre: The 1 hour observation with no errors. Right: A 1 hour image with 
gain/phase errors. In each case, the peak of the difference is approximately 
10% of the peak in the image, or ~ 20 mJy. 



cal RMS noise of 7.5/iJy beam ^.^ The best-fitting parameters for 
the lens model w as found using the downhill simplex method (e.g. 
IPress etalJl992l chapter 10). 

Before discussing the results of the simulations, there is an im- 
portant question that must be addressed: how do we know we have 
the best model? In other words, is there a well defined minimum in 
which is defined as the best model? To explore this question, an 
initial investigation was undertaken with a single simulated image. 
Parameter space was searched to understand the properties of the 
surface as a function of parameter value. 

Figs. ll3lfT4l and fT5l show slices through space for the criti- 
cal radius, ellipticity, orientation angle, core radius and lens centre. 
In each case, the surface varies smoothly with parameter value and 
has a well defined minimum. The best-fitting parameters are thus 
defined unambiguously by the minimum in the surface^. In ad- 
dition, the smooth nature of the surface and well defined minimum 
means we can be confident that the search algorithm will find the 
global minimum in x^- 

It is apparent from Figs. I13II14I and [T5l that there are some 
mild degeneracies between fitted parameters. There is a degener- 
acy between critical radius and core radius, and weak degeneracies 
between critical radius vs ellipticity and critical radius vs position 
angle. The degeneracy between critical radius and core radius is to 
be expected because a non-singular core can only be distinguished 
if it produces central images or changes the local slope of the mass 
profile enough to distort the inner image. In the case of the simu- 
lated lensed lobe, the bright inner image is 0.7 arcsec from the 
centre, hence it is reasonable that the models produce the same re- 
sult for small (< 0.1 arcsec) core radii. The x^ surface in Fig. ll4l 
slopes away to the right because the overall mass scale (i.e. b) must 
be increased to compensate for the mass 'lost' in the centre of the 
lens due to the non-singular core. However, the actual critical ra- 
dius and total mass within the critical radius remain the same. In 
the case of critical radius vs orientation angle, the effect is very 
small and probably due to the particular geometry of this simu- 
lated lens (i.e. it is essentially a double). If the lens was a quad or 
contained more bright regions it is difficult to see how the position 
angle could compensate for changes in the critical radius. 

For the case of critical radius vs ellipticity, again this is a small 
effect. A small change in critical radius can be compensated for by a 
change in ellipticity for two-image systems. This moves the critical 
line which can move images relative to the critical line. However, a 
ring image such as that of the simulated lobe contains a great deal 
of information about the ellipticity. Hence it is well constrained as 
seen in Fig |13l 

For the lens centre, the major axis of the contours are 
aligned with the angle between the centre of the source plane and 
the bright regions in the image. The effect on the images of moving 
the lens centre can be partially corrected by moving the source. 

Having shown that the surface has desirable properties, we 



using the same series of steps as described in the previous section, 
including self-calibration. 

The lens was modelled as a PIE? with all parameters free to 
change. Given the artifacts generated in the CLEANed images, in 
particular the apparent offset of the reconstructed image compared 
with the true image, the simulations were performed twice: once 
with the lens centre free to change and once with the lens cen- 
tre fixed. The PSF was taken directly from the Gaussian restor- 
ing beam and the noise in each pixel was set to the fixed theoreti- 



° It should be noted that the uncertainty in pixels in a CLEANed image is 
not well determined and the theoretical RMS noise is certainly a lower limit 
to the true noise in the pixels — especially in the brighter regions. Never- 
theless, without a solid theoretical understanding of the noise in CLEANed 
images only ad-hoc schemes can be used to model pixel noise in image 
space. Hence, using the theoretical noise is a good starting point. An al- 
ternate value is the off-source RMS noise in the CLEANed map, which is 
approximately 55/iJy beam~^ for the 1 hour simulations. 
^ One should always bear in mind the possibility of degeneracies between 
more than 2 parameters. In this case, a mass slope parameter is not being 
used so marginalisation over this parameter is not necessary 
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Figure 13. Cuts through the surface of the simulated lensed lobe. Parameters which do not vary are fixed to their best-fit value. Contours increase by 5%, 
10%, 15%, 20% and 25% from the minimum. Left: Ciitical radius (b) vs ellipticity (e). Right: critical radius (b) vs position angle (d). 
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Figure 14. Cuts through the surface. Contours increase by 5%, 10%, 15%, 20% and 25% from the minimum. Left: position angle (0) vs ellipticity (e). 
Right: Critical radius (b) vs core radius (rc). 





True 


Lens centre free 


Lens centre fixed 


Parameter 


value 


Mean 


std dev 


Mean 


std dev 


b (arcsec) 


1.35 


1.3446 


0.0008 


1.3511 


0.0004 


e 


0.065 


0.0644 


0.0011 


0.0668 


0.0004 


6 (degrees) 





-0.55 


0.32 


0.604 


0.14 


Tc (arcsec) 





0.0018 


0.0080 


0.0026 


0.0048 


xi^ (arcsec) 





0.0402 


0.0028 


fixed 




j/£ (arcsec) 





0.0070 


0.0026 


fixed 





Table 5. Mean and standard deviation of model parameters for 100 Monte- 
Carlo simulations based on the best-fitting PIEP model. 



now examine the results of the simulations. Table|5|shows the mean 
and standard deviation of best-fit parameters for the 100 simulated 
images. The results with the lens centre both free and fixed are in- 
cluded. The main result from the simulation is that the lens model 
has been recovered well but there are biases in some parameters, in 
the sense that the standard deviation of the fitted parameter value 
is less than the difference between the mean and true parameter 



value. This bias is a systematic error introduced by analysing the 
processed radio images rather then the measured visibilities. Fortu- 
nately, the relative bias is small: ~ 0.5% for the critical radius and 
~ 2% for the ellipticity. Compared with the error in the fitted pa- 
rameters (as described in the next section), these systematic errors 
are very similar in magnitude, and therefore do not make the results 
of the analysis unusable. It is encouraging that the lens systems can 
be modelled using the processed radio images with modest system- 
atic errors. The main properties of the system: critical radius, ellip- 
ticity and position angle have been recovered. The location of the 
lens centre suffers from the worst systematic error but the deviation 
is of the order 1 pixel width (0.05 arcsec in these images). 



3.2.3 Reconstruction significance and error estimates 

Having shown that the software can correctly reproduce the lens 
model parameters, this section focuses on quantifying the power of 
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Figure 15. More cuts through the x surface of the simulated lensed lobe. 



Model 



Best parameter values 



PIEP 20517 

SIS+7 36454 

SPEMD 85155 

SPEMD 27627 



1.352, e = 
1.358, 7 
1.446, q -- 
1.285, q -- 



■■ 0.0669, e ■- 
-- 0.0589, e 
-- 0.867, 9 = 
: 0.750, e = 



= 0.65 
= 0.77 

4.14, /3 = 0.8 
-2.44,^9 = 1.2 



Table 6. Best-fit model parameters and for a radio lobe simulation with 
alternate (incorrect) mass models. The SPEMD models had fixed /3 and 
7-c = 0. 



the software to distinguish between lens models, and the accuracy 
of the model parameter measurements. 

Here we test if data with the same assumed errors and default 
source pixel value (A) can pick the correct lens model from a num- 
ber of similar candidates. The data was tested using PIEP, SPEMD 
and SIS-1-7 models for a single noisy simulation. The simulations 
were performed with an assumed pixel noise of 7.5/^Jy and A = 
O.OlmJy. For the SPEMD, two values of the mass profile where 
used: P = 0.8 and /3 = 1.2. 

Table |6| shows the best-fit and lens model parameters for 
each model. The table shows, as expected, the incorrect mass mod- 
els generate worse values. The evaluation is qualitative but sup- 
ports the idea that radio images still contain much information 
about the lens when processed in image space. 

The next task in modelling the simulated lensed lobe is to 
quantify the uncertainty in the fitted parameters. If the noise in the 
image pixels is uncorrelated, then the uncertainty in each parameter 
can be calculated by examining the shape of the x'^ surface around 
the minimum. For radio maps, this is not immediately possible be- 
cause the noise in neighbouring image pixels is correlated by the 
beam, hence a change in the source which affects a single pixel in 
the image is spread over an entire beam. In this case the total flux in 
a beam is 22 pixel units, so a change which generates a Ax^ ~ 1 
in a given image pixel (with no beam) makes Ax^ ~ 22 over the 
real image (to first order). 

Using Table |5| the standard deviation of the parameter value 
is interpreted as the la error. With the known deviations and the 
well behaved properties of the surface around the minimum, 
one might ask whether the errors can indeed be derived directly 
from the surface. As in the previous section, confidence limits 
can be derived from the by marginalising over parameters. The 
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Figure 16. Marginalised contours of x^/beam area for the test lensed lobe 
for two parameters. Contours are 2.3, 4.6, 9.2 and 18.4 above the minimum. 



pixels are not independent in the radio map, so x^ should be di- 
vided by the integrated beam so that it reflects contributions from 
individual pixels in the map. Fig. ll6l shows the marginalised scaled 
(X /22) surface for the PIEP model with contours 2.3, 4.6, 9.2 and 
18.4 above the minimum, corresponding to confidence intervals of 
68, 90, 99 and 99.99 percent respectively. It is immediately appar- 
ent that the 68 percent interval corresponds very well with the la 
errors derived from the Monte-Carlo simulations. This was true for 
all parameters, indicating again that the properties of the x^ sur- 
face are suitable for error estimation. Hence by marginalising over 
parameters and dividing x^ by the total of the beam, the x^ surface 
provides a good estimate of la errors. This is extremely useful and 
means that lengthy Monte-Carlo simulations are not necessary. 



3.2.4 Source Reconstruction 

In Fig. ^1 the reconstructed lobe is shown. The left image is the 
original unlensed lobe. The centre image shows the reconstructed 
lobe from a 1 hour noiseless integration. The reconstructed lobe 
is very similar to the known true source even though the lensed 
image is a CLEANed map. The image on the right is the recon- 



© 2002 RAS, MNRAS 000.171021 



16 R. B. Wayth and R. L. Webster 




-10 12-1012-1012 
Arcsec West Arcsec West Arcsec West 



Figure 17. Reconstructions of the source for the simulated lensed lobe. Left: The original (unlensed) lobe. Centre: The idealised case where the visibilities 
contained no noise. The reconstruction is almost identical to the original with some low level structure around the brightest spot. Right: The reconstruction for 
the best-fitting PIEP model with real (noisy) data. In all images, contours double from 80 /.tjy beam~^ 



stmcted source for the best-fitting PIEP model from a 1 hour inte- 
gration with noise. Again, the lobe has been recovered quite well, 
with some low-level structure around the brightest spot. 

4 EXAMPLES 

In this section we use archival HST and VLA data to model the 
sources HST J15433-I-5352 and MG154 9-H3047. Lensview has 
also been used in a detailed analysis jWavth & W ebster 20031: 
IWavth et alJ l2005l) of the optical Einstein ring ER 0047-2808 
jWarren et alll996lll999l:lKoopmans & TreJ2003h . 

4.1 The optical lensed arc HST J15433-h53S2 

HST J15433-I-5352 i Ratnatunea etal."l999') was one of several 
gravitationally lensed objects found in the HST Medium Deep Sur- 
vey. The data are available (in several filters) as a WFPC2 asso- 
ciation from the HST archive including variance files. The pixel 
scale of the data is 0.1 arcsec pixel"^. In this analysis we use only 
the F450W image which has the best signal in the arcs and the 
faintest lens galaxy. A section of the F450W image containing the 
lens galaxy (L) and companion galaxy (G) is shown in Fig. ll8l 

The lens galaxy was subtracted from a 31 x 31 section of the 
data using the Galfit code (Pena et al. 2002) v2.0.3 using a de 
Vaucouleurs light profile. The lensed arcs were masked to prevent 
them from affecting the subtraction. The data, mask and galaxy 
subtracted image are shown in Fig. 1191 

Redshifts for the lens {z = 0.497) and source {z = 2.092) 
have been measured by Treu & Kooomans ( 2004). The lens galaxy 
has a nearby (4.7 arcsec, z — 0.506) companion and is likely part 
of a larger group so there will be a significant external shear in this 
system. The properties of the len sing galaxy and companio n galaxy 
have recently been measured bv lTreu & Koopmanj i2004h . In this 
analysis we will focus only on the properties of the overall mass 
distribution in the galaxy and how well it can be constrained using 
lensing. 

We define an annulus-shaped mask of 121 pixels to encom- 
pass the lensed arcs and surrounding few pixels of the image. Initial 
experiments found what seemed to be very good fits although the 
formal was poor. This prompted a comparison of the noise prop- 
erties of the image and the variance file. Using a 50 x 40 square 
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Figure 18. HST J15433+5352 lens galaxy (L) and companion galaxy (G) 
from the HST F450W image. 




Figure 19. HST J15433+5352 data. Left: The 31 X 31 cutout from the 
HST F450W data. Centre: The mask used when fitting the galaxy biightness 
profile. Right: the galaxy subtracted data. 



pixel section of blank sky near the lensed image, we generated a 
histogram of the data values. The blank sky had mean 0.06 and 
standard deviation 1.36 counts. The median a of the variance file 
over the same region was 0.754, indicating that the variance had 
been underestimated. To rectify this problem, 0.6 was added to the 
sigma image (i.e. the square-root of the variance) and a new vari- 
ance image was generated. 

To model the lens in HST J15433-I-5352, we again used a 
pseudo-isothermal elliptical potential, PIEP, with an external shear. 
The centre of the lens was fixed to the centre of the light profile 
and we used a zero core radius. The external shear is defined by its 
magnitude 7 and position angle 9^ . There are therefore 5 parame- 
ters (b, e, 7, 6-y) in the lens model. The companion galaxy, (G), 
and nearby group will contribute to the lensing properties of this 
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Model 



PIEP 



X params 

75.9 6 = 0.525 ± 0.015 
e = 0.16 ± 0.01 
= -120 ± 3 
PIEP + 7 73.1 b = 0.525 
e = 0.11 
0, = -134 
7 = 0.06 

= -104 

Table 7. Results for simple mass models of HST J15433+5352. For the 
PIEP + 7 model, there were many combinations of €,8^,^,9^ that yield 
very similar x^, but none that were a substantial improvement on the PIEP 
model. For this reason, no errors are given. 



system. We aim to determine if the data in an extended system such 
as this can distinguish betwe en lens models with and wi thout exter- 
nal shear. We note also that lTreu & Koopmansl i2004) performed 
an analysis of this system as part of a larger set of lenses using a 
completely different method. We compare results below. 

Notwithstanding the discussion about degrees of freedom in 
the source model in ^23\ the image plane mask of 121 pixels does 
not leave many degrees of freedom in the source before the tar- 
get becomes zero or negative. Initial experiments showed that 
source plane pixels had to be at most 1/2 the size of image pix- 
els and that the source is elongated. We chose the smallest possible 
source plane that could accommodate the reconstructed source with 
a source-to-image plane pixel scale ratio of 1/2. The final source 
plane is 10 x 10 pixels and offset from the centre of the lens plane. 
Using equation|8| the target is 18 ± 6. 



4.1.1 Results 

Best fitting parameters of the modelling are shown in TableQ The 
significance of the values is discussed in the next section. The 
surprising result from the modelling is that a purely elliptical model 
is equally able to reproduce the data as a model which includes 
shear. Put another way, the combined effect of the external pertur- 
bations plus intrinsic ellipticity in the lensing galaxy is indistin- 
guishable from a different, purely intrinsic, ellipticity (and orien- 
tation) in the galaxy itself (at least in this lens). This result is in 
contrast to other results (e.g. .Keeton et a l. 1997), which find that 
many QSO lenses with only positional information can only be 
explained with the presence of an external shear. A possible rea- 
son why an external shear is not required, at least in this lens, is 
that the lens is essentially a cusp image. The reconstructed source, 
shown in Fig. 1201 lies across a cusp. The lensing properties of 
cusps have been shown to be quite generic iSchneideret alll992t 
ISchneider&Weisslll992h so the section of the arc that lies very 
close to the critical line will be locally independent of the partic- 
ular details of the lens model. The section of arc that is not close 
to the critical line has its shape determined both by the lens model 
and the shape of the source, so it caimot reveal anything conclusive 
about the lens without prior knowledge of the source properties. It 
is reassuring, however, that a simple lens model and morphologi- 
cally uncomplicated source can explain the image. 

Errors on the fitted parameters are also shown in Table |7| for 
the PIEP model. The errors were generated using the properties of 
the x^ surface around the minimum. For the PIEP -I- 7 model, the 
ellipticity and orientation angle parameters e,9^,'y,9^ are virtu- 
ally degenerate hence no errors are given. Although our modelling 
found a 'best-fit' for this model, many combinations of parameters 
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Figure 20. HST J15433+5352 models. Left: the galaxy subtracted data. 
Centre left: The mask of 121 pixels used in the modelHng. Centre right: 
The best-fit PIEP model image and critical line. Right: The model source 
and caustic. 



produced x^ values which were within Ax^ ~ 3 of the best fit, 
therefore equally as good. For any (reasonable) set of values (e, 6'e) 
there exists (7, 9-y) such that the x^ is as good, but not significantly 
better, than the PIEP. 

The critical radius, 6, differs from lTreu & K oopmans^ i2004h 
(0.36 arcsec), however they modelled the companion galaxy, G, as 
an isothermal sphere. Given its close proximity to the lens galaxy, 
L, the halo of G will contribute substantial convergence to their lens 
model, so their smaller Einstein radius is merely making up the 
difference. The critical radius is also substantially different from 
the value found by Knudso n et al. (_2001) {b = 0.58 arcsec) who 
used simple elliptical isothermal models like those used here. This 
systematic difference is huge (10%!) compared to the random error 
from the fitted models. 

The model image and reconstructed source for the PIEP model 
is shown in Fig. |20| The PIEP -1- 7 model looks identical, so it is 
not shown. Our model shows a source, approximately 0.5 arcsec in 
length, which straddles a cusp. The brightest region in the source is 
outsi de the caustic. The source reconstruction is qualitatively simi- 
lar to lTreu & Koopmansl |2004) and Knudson et al. ( 2001) but dif- 
fers in detail. Comparing to lTreu & Koopmans, (.2004') . their source 
has the brightest region outsid e the caustic but the m odel is curved 
around the cusp. Compared to ' Knudson et alj i200lh . their source 
is placed at a different location, has a different position angle rel- 
ative to the caustic and is more round. The model image gener- 
ated by their source has significant residuals as seen in fig. 4 of 
lKnudsonet^ j200lh. In particular, the fainter part of the arc is too 
bright and the brightest region of the arc (at top and left of Fig. 1201 
LHS) are too faint. Given a parametric source model, however, it is 
the best that can be do ne. Hence, the use of a parametric source by 
lKnudsonet"^<200ll) has resulted in a significant systematic error 
in the critical radius of their best lens model. 



4.1.2 Source plane pixel thresholding 

At face value, neither of the models produce a statistically accept- 
able result. The reason for this is obviously the large source plane 
and the number of degrees of freedom used in the model. From Fig. 
I20l it is clear that the source pixel threshold method is appropriate 
for this model. 

Using this scheme, we must find T for our source. This is 
straightforward and we find setting T=2.5 makes 50 pixels go to 
zero and produces a Ax^ of 3.0 as shown in Table|8| The pixels 
set to zero are shown in Fig. 1211 Thus, our revised source plane 
has 50 degrees of freedom so the target x^ is 68 ± 12. Using this 
revised definition, the original models are acceptable. For compar- 
ison, the typical noise level in the image is 1.4 counts, the default 
source pixel value parameter. A, is 2.0 and the peak count in the re- 
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Threshold 


Number of 


Adjusted 




T 


zeroed pixels 


9 

target x 


New X 


1.0 


30 


48 it 10 


74.8 


1.5 


40 


58 ± 11 


75.5 


2.0 


46 


64 ± 11 


76.9 


2.5 


50 


68 ± 12 


78.9 


3.0 


55 


74 ± 12 


82.8 



Table 8. Results of adjusting the target and source for the PIEP model 
based on a threshold, T, below which source pixels are set to zero. In this 
example, the noise in each pixel is approximately 1 .4 counts and the default 
pixel value parameter, A. is 2.0. For each threshold level, the new source 
was re-projected through the best-fitting model and a new calculated. 
This table shows setting pixels with value < 2.5 to zero makes no difference 
to the model. 



Figure 21. HST J15433-I-5352 source models and source threshold masks. 



constructed source is 77.0. Hence we really are zeroing only those 
pixels which have not contributed anything to the model. 

Our method shows that the combination of using a mask on 
the image plus a source plane zeroing threshold leads to a sensible 
and achievable definition for the target ■ Formally, we are simply 
replacing Nsovltcc in |8j with Nsomcc > T. 

4.2 MG1549-I-3047 

MG1549-I-3047 (also known as 4C-I-30.29) is a. z = 1.17 radio 
lobe being lensed by a z = .11 SBO galaxy jLeharetaljl 19931 
1 19961 : IXreu & K oonmans 2003). We use previously unpubli shed, 
archiv al VLA data of MG1549-I-3047 in this analysis. Lehar et alj 
il993h used VLA data taken on 1991 August 2 in A configuration 
with a 20 min integration and achieved an RMS noise of 0.01 mjy 
beam~\ MG1549-^3047 was re-observed on 1992 November 18 
at 1.46 GHz (5 min), 4.89 GHz (30 min), 8.41 GHz (50 min) and 
14.96 GHz (190 min) again with the VLA A configuration. We use 
this new 8.41 GHz data in this work because the 14.96 GHz data 
has very poor signal in the ring. Three antennas (10,27,29) were 
off-line during the observations. 

We followed a standard VLA calibration procedure using 
AIPS and calculated 3C286 fluxes using the AIPS 'setjy' task. The 
data were quite clean and we used one round of self-calibration 
to make an image. The beam size was 0.21 x 0.24 arcsec. The 
off-source RMS noise in the map is 0.041 mJy beam~^ which is 
almost twice the theoretical no ise of 0.025 mJy beam^ ^ and four 
times the RMS noise stated in iLehar et alj 11993) table 2. Since 
fig. 2 in iLeharetalJ 11991 shows contours doubling from 0.2 mJy 
beam~^, there may have been a typographical error and the ac- 
tual noise in their map was 0.13 mJy beam~^ rather than the stated 
value of 0.013 mJy beam~^. The larger value would also be much 
more consistent with the theoretical value of 0.05 mJy beam^^ for 
a 20 min integration. 

HST NICMOS F160W and WFPC2 image associations (pro- 
posal id: 7495, PI: Falco) were obtained from the HST archive 
which show the faint z = 1.17 radio galaxy, which is the source 
of the lobes plus a nearby z = 0.604 galaxy iTreu & Koopmanl 



Parameter 


This 
work 




L93 
result 


b (arcsec) 


1.105 


±0.005 


1.15 




0.056 


±0.006 


0.073 


0e (deg E of N) 


-47 


± 1 


-48 


s,Tc (arcsec) 


0.132 


±0.03 


0.22 


XL (arcsec E of radio core) 


-6.59 


±0.02 


-6.64 


yL (arcsec N of radio core) 


3.33 


±0.02 


3.38 


Mass enclosed (x 10^" M©): 


7.06 


±0.06 


7.72 


cr„ (for best fit b) (km s~^) 


209 




213 



Table 9. Best fit PI EP model paramete rs for MGI549-^3047 and coiTe- 
sponding values from lLehar et alj il993h . converted to our cosmology. 

l2003t) . We used the HST images purely for astrometric purposes 
and performed no processing other than that which is automated by 
the archive pipeline (flat fields, CR rejection, image combining). 
We aligned the radio and optical images using the radio galaxy 
with a combined positional uncertainty of 0.04 arcsec. The com- 
bined radio and optical images are shown in Fig. 1221 We calculate 
the centre of the lens galaxy to lie at (J2000) RA: 15''49'"12^330, 
DEC: 30°47'16'.'51 with an uncertainty of 0.04 arcsec. The ring 
has lensed flux almost all the way to the centre of the galaxy, so it 
is extremely important that the location of the lens centre is known 
accurately. If the location of the lens centre is incorrect, any mea- 
surement of the properties of the mass in centre of the galaxy will 
be affected. 

We model MG 1549-1-3047 using a simple PIEP lens model to 
com pare to previ ous work and derive some basic properties of the 
lens. Leha r et"ai] jl993 !h presented a lens model based on matching 
various features of the lensed image, then performing a QSO-like 
analysis using the locations of the features. Their lensing potential, 

= (3(jl + {l-e){x/sf + {l + e){y/sfY^^-l^ has the 

same deflection angle properties as the PIEP with b = f3/s and 
rc = s. 

W e searched parameter space using the results from 
iLehar e t al. 1 1993) as a starting point. The parameters were unre- 
stricted except for the lens centre, which we allowed to vary by up 
to 0. 1 arcsec from our measured position. A 99 x 99 grid of pixels 
from the radio data was used for the modelling, shown in Fig. 1231 
(LHS). The variance in each pixel was defined to be (0.04)^, i.e. 
the square of the off-source RMS noise in the map. 

Th e result of our mod elling, along with some results from ta- 
ble 4 in lLeharetalJiT993) are shown in Table |9| The errors were 
calculated as discussed in il3.2.3l The most significant difference is 
that of the Einstein radius. Since the mass enclosed within the ring 
is oc b^, the difference in Einstein radius changes the mass by 10 
percent, which is quite a substantial amount given the accuracy to 
which Einstein radii are usually measured. Our predicted velocity 
dispersion, 209 km s~^, is slightly lower, but consistent with, the 
measured value of 227 ±18 km s" fLeharetalJl993) . 

The lens centre was found to be very close (0.024 arcsec) 
to the position predicted by the combined optical/radio astrome- 
try. The best fitting lens centre is at RA: IS** 49™ 12330, DEC: 
30°47'16'.'42. The uncertainty on the fitted lens centre (0.02 arcsec 
in each direction) is smaller than the uncertainty from the combined 
astrometry. This result shows that the software can be used to find 
the centre of the lens in cases where there is some uncertainty. We 
found that the best-fitting model core radius, s, was affected by the 
position of the lens centre. This can be understood by noting that 
as the lens centre moves East of its known position, a radial critical 
line is no longer required to truncate the image. Rather, the image 
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Figure 22. MG1549+3047 radio galaxy and lobes. Left: 8.4 GHz radio contours are overlaid on the HST WFPC2 F814W image. Right: 8.4 GHz radio contours 
are overlaid on the HST NICMOS F160W image. Radio contours double from 0.2mjy. Between the radio core and lens is an intermediate redshift (z=0.604) 
galaxy. 
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Figure 23. MG1549+3047 data and models. North is up. East is left and position (0,0) represents the centre of the lens. Left: The lensed lobe. Centre: Model 
PIEP image with tangential critical line (dashed). Contours increase by from 0.2mjy beam^^ for the data and model. Right: The reconstructed source 
after smoothing with a boxcar filter of width 3. The caustic is also shown. Contours for the source increase by 2 from 0.08 mjy pixel~^ 



is simply disappearing at the centre of the lens. Despite the cor- 
relation between parameters, the lens centre was unambiguously 
identified by the modelling. 

Another noteworthy property of the model is the requirement 
of a small constant surface-density core. The lensing galaxy has 
a prominent bulge which should follow an r^^^'-type light profile. 
(We leave a detailed analysis of the properties of the lens galaxy 
to a future work.) Such a profile has well defined lensing deflec- 
tion properties. The magnitude of the deflection rises from zero in 
the centre of the galaxy, reaches a peak around 0.2 scale radii then 
falls away roughly as 1/r, as shown in Fig. 1241 Any bulge-like 
galaxy's gravitational and lensing potential should be dominated 
by stars in the centre of the galaxy, so one would expect the lensed 
image to reflect the lensing properties of the bulge in the central re- 
gions of the galaxy. Further out, the galaxy's dark matter becomes 
more and more important so the lensing properties follow the fa- 
miliar isothermal model. Our model confirms this expectation with 



the requirement of a core (non-isothermal) region with zero/small 
deflection angles. 

The reconstructed source, shown in Fig. 1231 looks just like 
a regular FR-I type radio lobe. It has many similarities with its 
unlensed counterpart, including a hot-spot, diffuse tail and some 
slightly brighter structures in the tail. The total flux in the recon- 
structed lobe is 35mJy compared to 44mJy in the unlensed lobe. 
The overall magnification for the lobe is modest: about 3.4, but this 
is because the bright spot in the lobe is outside the multiple imag- 
ing region of the lens. It is the much fainter tail which is the most 
highly magnified. 

The significant difference in Einstein radius between our 
model and L93 highlights the importance of using detailed mod- 
els for resolved lenses. The analysis of L93, based on identifying 
features in the lens, is dangerous- especially when identifying fea- 
tures close to the centre of the lens in low magnification regions. 
L93 identified three sets of features (AI,A2,BI,B2,CI,C2) in the 
lens. Our model shows that A2 and B2 are features which have 
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Figure 24. Relative magnitude of the deflection angle for a de Vaucouleurs 
model with the radius in units of the galaxy 1/2 light radius, R\/2- 

several beam areas worth of source brightness information packed 
into one or two beam areas in the image. Since the source bright- 
ness varies rapidly with position in the lobe, the resulting brightness 
at position A2 and B2 cannot be directly mapped to another region 
of equal brightness in the image. Al and Bl are a good approxima- 
tion of the true surface brightness in a particular part of the source, 
whereas A2 and B2 are the average brightnesses of a large portion 
of the source spanning several beam areas. 



5 DISCUSSION AND CONCLUSIONS 

Resolved gravitationally lensed images provide the power to dis- 
criminate between mass models for the lensing galaxy. However, 
until recently this potential has been largely untapped. There is now 
a growing set of algorithms and tools for modelling resolved lensed 
images. In this section we discuss issues relating to our software, 
Lensview and how it compares to other algorithms. 

Mapping matrix vs non mapping matrix: The mapping matrix 
has many convenient properties, not the least of which is that it 
correctly preserves surface brightness when projecting between the 
source and image planes. Warren & Dve (2003) use ray-shooting 
to map source pixels to image pixels. This leads to the question: if 
a projected image pixel overlaps with several source pixels, which 
one do you choose? This is most important for image pixels which 
are in a medium to low magnification region. Such pixels will over- 
lap with several source pixels and their true value has a contribu- 
tion from all. By ray shooting, only a single pixel is sampled. The 
effects of this sampling are unknown at prese nt as are the effects 
of the simple interpolation scheme proposed bv lTreu & KoopmanJ 
i2004h . In principle, there is no reason why a mapping matrix 
caimot be incorporated into the semilinear method, but the in- 
creased source-to-image pixel relationships of the mapping matrix 
means the matrix be to inverted by the semilinear method rapidly 
loses its sparseness. Obviously a mapping matrix is not useful for 
LensCLEAN which uses point sources, although the LenTil al- 
gorithm I Wucknitz 2004) uses something similar to find image po- 
sitions. 

Iterative vs non iterative: The source model for our method 



builds up through progressive iterations of forward and reverse pro- 
jection. The semilinear method calculates the source model in a sin- 
gle matrix inversion where the matrix has size N x N (with A*' pix- 
els in the image) and is presumably fairly sparse. The 'sparseness' 
of the matrix must be reduced with an increasingly broad PSF. This 
will affect running times and possibly be susceptible to accumu- 
lated numerical errors. The iterative approach is not significantly 
affected by the size of the PSF and numerical errors do not ac- 
cumulate between iterations. The time to project between source 
and image planes, however, is affected by the fraction of the image 
plane which covered by the source plane. Nevertheless, the major- 
ity of the computational time in our method is taken convolving 
with the PSF. Typically 100-200 iterations are required. 

Running times: Code execution times have not been given for 
other methods. As mentioned above, LENSVIEW spends most of its 
time performing convolutions via the FFT Hence running times are 
roughly 0{N log N) where A'^ is the number of pixels in the im- 
age. For the examples in this paper, running on a 2.0GHz Pentium 
computer, the inner loop is completed in less than one second for 
J15433 and in approximately seven seconds for MG1549. Hence, 
large sections of parameter space can be searched in reasonable 
times with LENSVIEW. 

Entropy/regularisation: The entropy constraint is analogous 
to regularisation of semilinear inversion although it should be em- 
phasised that the way in which the source is reconstructed in each 
method is quite different. Any form of regularisation introduces 
bias into the source, however that is the price one must pay for us- 
ing real (noisy) images since all image deconvolution methods will 
require regularisation in some form. The entropy constraint is an 
excellent form of regularisation because it smooths source pixels 
which have little or no data constraining them and enforces positiv- 
ity in the source. Other choices for a source constraint are possible, 
but we have not explored them here. 

Images vs visibilities: Although LENSVIEW currently works 
only with images, the extension to use radio visibilities is conceptu- 
ally straightforward. It is difficult to see how the semilinear method 
could be extended to use radio visibilities because the matrix inver- 
sion would become impractically large- this is exactly why CLEAN 
and similar algorithms are used to process radio images. 

A simple estimate shows that processing the data in image 
space provides approximately two orders of magnitude increase in 
the speed of analysis. Fourier transforming between the visibilities 
and dirty images used in the simulations takes ^ 1 second. Iterative 
methods must be used with this volume of data, and our models 
(and those of SB84) usually converge in 100-200 iterations. As- 
suming no other computational burden, the forward and backward 
Fourier transforms alone will take several minutes for a single lens 
model for the inner loop. Conversely, by using the CLEANed im- 
age, just the region of interest can be extracted from the larger im- 
age for use with LENSVIEW. 

The tests presented here show that the accuracy of the lens 
model parameters is biased by working in image space, however 
the systematic errors are sm all enough that t he res ults are still 
useful for some applications. lEllithorpe et alj il996l> found sim- 
ilar biases after comparing results of image-space models with 
visibility-space models. They the refore concluded that using vis- 
ibilities is preferable to images. IWucknit3 i2004b subsequently 
show ed that further improve ments to LensCLEAN were possible. 
Since lEllithorpe etlll jl99d t) did not perform any simulations, it is 
not known whether their visibility-based result s were actually af- 
fected by the problems identified by Wuck nitJ i2004h . Hence it is 
possible that systematic effects were still present. 



© 2002 RAS, MNRAS QOQ.ITEl 



Lens VIEW 21 



A major advantage of using visibilities is that a well de- 
fined target ex ists. In some rare cases, e.g. MGl 131+0456 
teewitt et all 1983) . the entire image generated by the interferome- 
ter is contained in a very small area — just a few square arcseconds. 
In that case it would be possible to process the dirty image directly 
using the dirty beam as a PSF. Since there would be no artifacts due 
to image processing in that case, the noise properties and target 
of the image would be well defined. 

Mass distributions: Lensview can work with an arbitrary 
lens mass model with negligible extra computational cost. This 
appears to be true also for the semilinear method, but is not so 
for LensCLEAN which requires solutions to the lens equation 
for each point source component. Hence LensCLE AN is substan- 
tially slower when using a non-analytic lens model. 

The potential applications of LENSVIEW extend well beyond 
the simple examples shown in this paper. Some of the most exciting 
applications are: 1) determining the slope of the mass profile in lens 
galaxies. Interpretation of Ho measurements rely on this knowl- 
edge. Wavth et al. ( 2005) have shown that the optical Einstein ring 
ER 0047-2808 has a well determined mass profile in the region of 
the images. 2) properties of the dark haloes of lens galaxies. It is 
straightforward to model a lens as a stellar (baryonic) component 
with a dark matter halo. Extended images will allow measurements 
of the properties of the halo (e.g. Wavth & Webster 2003). 3) deter- 
mination of the lens centre. Again, required for Ho measurements. 

We have presented simple elliptical isothermal lens models for 
the optical lensed arc HST 115433+5352 and the radio Einstein ring 
MG1549+3047. Our models show that for HST 115433+5352, the 
data are unable to distinguish between an elliptical model and an 
elliptical model with external shear. We speculate that this because 
the source is lying over a cusp, therefore has a very generic image 
configuration. For MG1549+3047, we show that the previous lens 
model had an incorrect Einstein radius and therefore the mass en- 
closed in the ring was overestimated. In addition, we showed that 
the model requires a small non-singular core, consistent with ex- 
pectations for a star-dominated mass distribution in the centre of 
the lens galaxy. 

More detailed applications of the software are possible. We 
look forward to a new era of precision modelling of resolved im- 
ages. 
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